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Abstract. This paper is the fourth one of a series whose 
chief objective is studying the influence of different mass 
spectra on the dynamical evolution of open star clusters. 
Results from several A-body calculations with primordial 
binaries and mass loss due to stellar evolution are pre¬ 
sented. The models show significant differences with those 
for primordial binaries but no stellar mass loss presented 
in de la Fuente Marcos (1996b). A differential dynami¬ 
cal behaviour depending on cluster richness is found com¬ 
pared to de la Fuente Marcos (1996a). The evolution of 
these realistic models is very dependent on the initial mass 
function. Even for rich clusters, there is a dependence on 
the binary mass spectrum. The velocity distribution of 
the escapers is examined and compared with results from 
previous calculations. The evolution of the primordial bi¬ 
nary population is analyzed in detail. The cluster remnant 
and the final binary population are also studied. Finally, 
some conclusions about observational properties of Open 
Cluster Remnants are presented. 

Key words: chaotic phenomena - celestial mechanics, 
stellar dynamics - Galaxy: open clusters and associations: 
general - stars: binaries: general - stars: evolution - stars: 
luminosity function, mass function 


1. Introduction 

In the last few years, a significant number of pre-main- 
sequence binaries and multiple systems have been discov¬ 
ered in young open clusters (Brandner et al. 1996; Ghez 
et al. 1993, 1994; Leinert et al. 1991, 1993; Padgett et al. 
1996; Prosser et al. 1994; Richichi et al. 1994; Simon et 
al. 1992, 1993, 1995). All the surveys carried out find a 
binary frequency which is greater than or equal to that of 
the field stars in the range of separations to which they are 
sensitive. The observational data suggest that the binary 


cluster population can be interpreted in terms of differ¬ 
ent formation mechanisms. Wide binaries (with periods 
greater than 100 years) could originate in capture events 
or by a fragmentation process during the collapse of a sin¬ 
gle rotating protostar. In fact, Hartigan et al. (1994) have 
shown that about one third of wide pre-main-sequence bi¬ 
nary pairs in their sample (projected separations of 400- 
6000 AU) are not coeval, with the less massive star usu¬ 
ally being younger than the more massive star, suggest¬ 
ing that they are formed by capture or exchange (more 
likely) events. On the other hand, close systems must al¬ 
most certainly be primordial, although their exact origin 
is not yet well understood. Several mechanisms have been 
suggested (fragmentation during the late collapse, gravi¬ 
tational instabilities or even orbital decay) for explaining 
the formation of these close primordial binaries. Moreover, 
observations indicate that binary stars are seen since the 
earliest stages of star formation, which suggests primor¬ 
dial mechanisms for their origin (Harjunpaa et al. 1991; 
Reipurth & Zinnecker 1993; Mathieu 1992, 1994, 1996). 
In any case, observational results hint that star formation 
produces mostly binaries. Further information about bina¬ 
ries in clusters can be found in Trimble (1980), Abt (1983), 
Reipurth (1988), Mathieu (1989), Zinnecker (1989) and an 
extensive recent review in Bodenheimer et al. (1993). 

Binaries in star clusters are of chief importance both in 
observational and theoretical astrophysics. In some open 
clusters, a certain number of stars appear above the cluster 
main-sequence turn-off (blue stragglers); it is currently in¬ 
terpreted (Wheeler 1979) as a result of stellar coalescence 
or extended main sequence life-times caused by mixing 
within the stars (Abt 1985). Another popular explanation 
for the origin of blue stragglers may be due to binary mass 
transfer (McCrea 1964). Also, runaway OB stars are inter¬ 
preted as binaries escaping from star clusters (De Cuyper 
1982; Sutantyo 1982; Hills 1983; Gies & Bolton 1986). In 
addition, binaries in open clusters are used to determine 
the distance scale (standard candles). Binaries are also the 
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key to explaining X-ray emissions from globular and old 
open star clusters. From a theoretical point of view, it is 
clear that primordial binaries (hereafter PBs) can domi¬ 
nate the early stages of the dynamical evolution of star 
clusters (see Heggie & Aarsetli 1992). 

Although binaries are heavier than the mean mass in 
a cluster, their masses do not remain constant throughout 
their life. Mass loss from stellar evolution, stellar collisions 
and exchange of mass from one component to another can 
alter significantly the binary orbit or even disrupt it. Pri¬ 
mordial binaries which are too close will be modified or 
even disrupted by mechanisms such as mass loss or mass 
transfer during the evolution of the primary from a main 
sequence star to a giant. Because of their massive nature, 
binaries tend to occupy the inner regions of star clusters 
so a change in their physical parameters (mass is the most 
important) may affect significantly the entire dynamical 
evolution of the system. Their impact on the cluster evo¬ 
lution depends for the most part on the coupling between 
the characteristic time-scale for mass segregation in the 
cluster and that for stellar mass loss. If mass segregation 
can be reached before significant mass loss has occurred, 
these primordial binaries act as energy sources more or 
less continually until they become unstable by stellar evo¬ 
lution or escape. In the case of important stellar mass loss 
before equipartition, this situation may not be attained 
and the new evolved binaries may have no tendency to 
reach the inner regions of the cluster because their masses 
have decreased in a significant way. However, from a the¬ 
oretical point of view, a smaller effect is expected than 
in the case of models without PBs but stellar evolution, 
because the single stars are also losing mass, so binaries 
are in any case the most massive objects. 

The role of stellar evolution for the dynamics of star 
clusters has been of interest since the end of the seventies 
(Angeletti & Giannone 1979, 1980; Angeletti et al. 1980) 
and there are a number of recent papers (Stodolkiewicz 
1982, 1985; Terlevich 1983, 1985, 1987; Applegate 1986; 
Weinberg & Chernoff 1988; Quinlan & Shapiro 1990; de la 
Fuente Marcos 1993, 1996a (hereafter Paper II); Aarseth 
1996b, c). 

The effect of PBs on the dynamical evolution of clus¬ 
ters was studied before stellar evolution. The first paper 
on the subject was by Aarseth (1975) and it was followed 
by a number of papers (Aarseth 1980; Spitzer & Math- 
ieu 1980; Giannone & Molteni 1985; McMillan et al. 1990, 
1991a, b; Murphy et al. 1990; Gao et al. 1991; Hut et al. 
1992; Heggie & Aarseth 1992; McMillan 1993; McMillan 
& Hut 1994; de la Fuente Marcos 1996b (hereafter Pa¬ 
per III)). In addition, Hills (1975) gave a semi-analytical 
discussion about this question, Goodman and Hut (1989) 
pointed out the importance of PBs for the evolution of 
globular clusters and Leonard and Duncan (1988, 1990) 
carried out ./V-body simulations with primordial binaries, 
although their main emphasis was on escapers. Summaries 
of most of these studies can be found in the introductory 


section of McMillan et al. (1990), Gao et al. (1991) and 
Heggie and Aarseth (1992). 

On the contrary, the study of the interplay between 
stellar mass loss and PBs on the evolution of star clus¬ 
ters has only very recently become of interest. Pols and 
Marinus (1994) studied the binary stellar evolution in 
young open clusters using Monte-Carlo simulations, al¬ 
though their chief interest was in pure stellar evolution, 
and not in the dynamical one. Direct ./V-body calculations 
have been performed principally by Aarseth, but only a 
few details have been published (Aarseth 1996b, c). He 
has found that because of stellar evolution the fraction of 
binaries increases in the central regions of rich star clus¬ 
ters (N up to 10 4 ). This increase in the central binary 
fraction is because massive single stars evolve to low-mass 
stars. 

This paper is mainly devoted to study the interplay be¬ 
tween the mass spectrum, the PB fraction and the mass 
loss due to stellar evolution on the dynamical evolution of 
open star clusters. Moreover, we want to compare our re¬ 
sults for the surviving binary fraction with observational 
data for binaries in open clusters. We are mainly con¬ 
cerned with the binary type for trying to answer the ques¬ 
tion about the preferential type of surviving binaries in 
open clusters. It is to be expected that binaries with both 
components being low-mass stars (late spectral types) will 
be preferential survivors in rich open star clusters because 
the time-scale for cluster disruption is larger than their 
characteristic time-scale for significant mass loss due to 
stellar evolution. However, for poorly populated open clus¬ 
ters, the disruption time-scale can be significantly smaller 
than the stellar evolution time even for moderately mas¬ 
sive stars, so it should be possible to find binaries with 
massive components. Also, we are interested in compar¬ 
ing the present results with those from previous papers in 
this series (de la Fuente Marcos 1995 (hereafter Paper I); 
Paper II; Paper III) concerning the role of the initial mass 
function (hereafter IMF) on the dynamical evolution of 
cluster models. 

We have performed five runs each for a total number of 
stars N = 100, 250, 500, 750 with five different IMFs using 
direct integration methods. As in previous papers, we use 
the same version of the Aarseth’s code NBODY5 (Aarseth 
1985; Aarseth 1996a). This code has become a standard 
in the field of star clusters simulations. Written in FOR¬ 
TRAN, it consists of a fourth-order predictor-corrector 
integration scheme with individual time steps. It utilizes 
an Ahmad-Cohen (1973) neighbour scheme to facilitate 
calculation of the gravitational forces, and handles close 
encounters via two-, three-, four-, and chain regularization 
techniques (Kustaanheimo & Stiefel 1965; Aarseth & Zare 
1974; Mikkola 1985; Mikkola & Aarseth 1993). 

All the calculations have been performed on a VAX 
9000/210, running under Open VMS operating system, at 
the Centro de Proceso de Datos (UCM, Moncloa, Madrid). 
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This machine has one CPU and its peak performance is 
about 100 Mflops. 

2. Cluster models 

In the present section all the physical features of the mod¬ 
els we have performed are discussed. This description fol¬ 
lows closely that in the previous papers of this series, 
although in this case all the astrophysical processes dis¬ 
cussed above are included simultaneously. 


Table 1. IMFs used in the calculations 

Salpeter IMF a = 2.35 

Taff IMF a = 2.5 (N < 100) 

a = 2.65 (N > 100) 

Miller & Scalo IMF (Eggleton et al. 1989)* 

Kroupa IMF (Kroupa et al. 1993) 

Scalo IMF (Eggleton 1994) + 

* This IMF is a fit to the Miller & Scalo’s (1979) results. 
+ This IMF is a fit to the Scalo’s (1986) results. 


2.1. Initial Mass Function 

The frequency distribution of stellar masses at birth is a 
fundamental parameter for studying the evolution of star 
clusters. As in previous papers (Paper I; Paper II; Paper 
III) several IMFs have been used in the calculations in 
order to generate an initial distribution of masses. We refer 
to Paper I for a full discussion of the IMFs used. Models 
for Kroupa (Kroupa et al. 1993) and Scalo (Scalo 1986; 
Eggleton 1994) IMFs are used with the binary correlation 
described in Paper III (Eggleton 1995). In Fig. [[] we see 
the modified Scalo IMF with the binary correlation used. 
For models with Salpeter (1955), Taff (1974) and Miller & 


Fig. 1. Distribution function for the modified Scalo IMF with 
the binary correlation described in Paper III. Note the loga¬ 
rithmic scale for both axes. 


Scalo (1979) IMFs the two components of the binary have 
the same mass. The five IMFs used in our calculations are 
summarized in Table [I]. 

2.2. Stellar evolution 

As in Paper II we use the fitting functions by Eggleton 
et al. (1989) in order to obtain the stellar diminution of 
mass as a function of time for Population I stars. These 
interpolation formulae are explained in the original Eggle¬ 
ton et al’s paper and partially in Paper II. All the stars 
start on a zero-age main sequence (hereafter ZAMS) with 
a uniform composition of hydrogen, X = 0.7, helium, Y 


= 0.28, and metallicity, Z = 0.02. For computational con¬ 
venience, the mass loss is implemented at discrete inter¬ 
vals (accumulated mass diminution of 1 %). Because of 
very different relative velocities, the actual mass loss is 
assumed to be instantaneous. The expelled gas leaves the 
cluster without any effect on the cluster members. We do 
not consider mass transfer processes in binaries, so the 
evolutionary times given by Eggleton’s fitting formulae 
are not changed. Also, there is no presence of accretion 
disks around stars, so disk accretion can not affect the 
evolutionary tracks of the cluster stars. There is only one 
way of changing the evolutionary time-scales for a given 
star; it can be achieved if two stars collide forming a new 
object, blue-straggler or Thorne-Zytkow object (Thorne 
& Zytkow 1977). 

2.3. Binary fraction and other binary parameters 

The number of primordial binaries in the cluster is conve¬ 
niently parameterized by the binary fraction / given by: 


/ 


N b 

N b + N s ’ 


(1) 


where N b , N s are, respectively, the number of stars which 
are binaries and singles. As in Paper III, a binary frac¬ 
tion of / = 1/3 is used in the calculations; so the overall 
multiplicity, defined as the ratio of the number of multi¬ 
ples to the total number of systems, is 0.33 in our present 
models. Lower bounds from observational surveys in Star 
Formation Regions are 0.37 (Ophiuchus) and 0.55 (Tau¬ 
rus) (Simon et al. 1995). For a sample of stars in the solar 
neighbourhood, Duquennoy and Mayor (1991) have found 
a fraction of 0.57, after correcting for observational bias. 
For comparing directly with models from previous papers, 
the number of singles and binaries is chosen in such a way 
that the total number of stars (not objects) is 100, 250, 
500 and 750 respectively. This arbitrary choice has minor 
effects on the cluster dynamical evolution as described in 
Paper III. 

In order to include a realistic initial binary popula¬ 
tion there are some other parameters to be determined 
in addition to the binary fraction. Our initial population 
of binaries is hard because they are of main dynamical 
importance for the evolution of star clusters. A binary is 
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defined to be hard or soft (Heggie 1980) depending on 
whether its binding energy is greater than or less than the 
local mean kinetic energy per star. As in Paper III, the 
semi-major axis of the binaries is taken from a uniform 
distribution: 

a b = a° b 10" 9 , (2) 

where a b is an input parameter whose value is about 1 /N 
in units of R v n (for a hard binary) and q is equal to X 
log72. The virial radius, R V i r , is defined by 

R vir = - GM 2 /4E , (3) 

where E is the total energy of the system, excluding the 
binding energies of any initial binaries, G is the gravita¬ 
tional constant and M is the total mass of the cluster; 72 
is an input parameter, and X is a random number uni¬ 
formly distributed in the interval [0,1]. The spread in the 
energy of the binaries is given by the spread in semi-major 
axis and finally it is given by 72. Small values of 7 Z pro¬ 
duce wide binaries, large values give close binaries. The 
value of 7Z used in our simulations is the same that we 
used in Paper III, 72. = 10. A typical value for the initial 
semi-major axis of the binaries in our simulations is about 
500 AU (the smaller value is about 70 AU and the max¬ 
imum is about 3300 AU). The values for the semi-major 
axis in Table || are the upper cut-off for the semi-major 
axis distribution. Eccentricities are chosen from a random 
(thermalized) distribution (Jeans 1929) and the same is 
done for the pericentre, node and inclination. The mass 
ratio for the PBs in all the models is 1/2, excepting those 
in which the binary correlation has been used. 

2.f. Main features of the models 

All the models (for the same N) have the same sequence of 
random numbers for generating initial conditions; spheri¬ 
cal symmetry and constant density are assumed, with the 
ratio of total kinetic and potential energy fixed at 0.25. 
Another characteristics common to all models are ran¬ 
dom and isotropic initial velocities. Also, in all the cases 
the cluster suffers mass loss due to escape of stars. A 
star escapes and is removed from the calculation when 
its distance from the cluster centre is greater than twice 
the tidal radius. The tidal radius is given by the clas¬ 
sical expression r t « (^fp) 1 ^ 3 where Ti is defined as a 
function of A and B Oort constants of galactic rotation: 
Ti = 2w[(— A — B) — w], where to is the rotational velocity 
of the star cluster around the galactic nucleus. As in pre¬ 
vious papers, the galactic gravitational field is introduced 
in the models as described in Terlevich (1987). Moreover, 
we ignore the effect of field stars on the dynamical evolu¬ 
tion of the cluster due to their high relative velocities (see 
discussion in Paper II). 

We use a standard and consistent set of units through¬ 
out, except where explicitly noted. All lengths are mea¬ 


sured in parsecs. Times are measured in terms of the ini¬ 
tial half-mass crossing time, defined by 

T cr = GM 5/2 /{-2E) 3/2 . (4) 

It represents the time taken for a typical star, moving with 
velocity < v 2 > 1 / 2 , to cross the virial diameter (2 /?„>). 
As in previous papers of this series we adopt a consis¬ 
tent tidal field for the models which have the same mean 
stellar density. A typical value for the mean mass density 
in a real cluster is 1.3 M@pc~ 3 (Lohmann 1971, 1976a, 
1976b, 1977a, 1977b) with a range of 0.5-3.2 Af@pc -3 
for his sample of open clusters. These values are typ¬ 
ical for evolved clusters, so that a larger mean density 
of ~ 12 Mp,pc~ 3 is adopted for the whole cluster (Table 
||) as considered young and not evolved. All the models 
do not have the same maximum (-M max ) and minimum 
(A! m in) masses for generating the IMF. The Salpeter,Taff 
and Scalo IMFs have A4 m ax = 15.0 Mq and _M m i n = 0.1 
M.q but the Kroupa and Eggleton IMFs use an algorithm 
that changes upper and lower limits for masses (see Table 
[|). Hence models for the Salpeter and Taff IMFs have the 
value of mean density quoted above but the others have 
~ 6 MqPc~ 3 (if using the same mean stellar mass, this 
will give the initial virial radius for different N). 

Table gives the disruption time for our present mod¬ 
els. If we compare these values with those from Paper III 
we observe that for N = 100 the disruption time is now 
smaller only for models with Salpeter, Taff and Miller & 
Scalo IMFs. However, for models with Kroupa or Scalo 
IMFs the disruption time is increased even comparing with 
Paper II. The main difference between the two groups of 
IMFs is the maximum mass. For the first group, the up¬ 
per cut-off for the mass spectrum is 15 solar masses but 
this value is significantly smaller for the second group, so 
the reason for this behaviour is clearly due to the mas¬ 
sive stars. It seems that the supernova events in the core 
destabilize the cluster in a very efficient way for poorly 
populated clusters. Two supernova events are enough to 
provide the energy source which disrupts the cluster. The 
acceleration of disruption as compared with Paper II is 
very significant for models with power-law IMF, however 
for Miller & Scalo IMF the disruption time in the current 
models is greater than the respective ones for Paper II 
(due to an upper cut-off in the mass spectrum). For mod¬ 
els with N = 250, the same trend pointed out above is 
observed for power-law models with regard to Paper III. 

On the contrary, realistic IMFs show greater disrup¬ 
tion times than those of Paper II and III. For N = 500, 
the behaviour of the disruption time is very similar to N 
= 250. For N = 750, power-law models have disruption 
times greater than those of Paper III but a bit smaller than 
those of Paper II. For Miller & Scalo and Kroupa IMFs the 
disruption times are smaller than those of Paper III and 
Paper II, however for the Scalo IMF the disruption time is 
greater than those from Paper II and III. This suggests a 
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Table 2. Main characteristics of the models 


MODEL 

IMF * 

N* 

N b 

a 0 x 

0 

■ML* 

M f ■ 

min 

< M >t 

Rl 

r t 0 

< R >o 

7■»© 

1 d 

I 

SA 

100 

25 

0.0160 

15.0 

0.1 



1.0 



1.26 

5.96 

0.78 

30.4 (74) 

II 

TA 

100 

25 

0.0160 

15.0 

0.1 



1.0 



1.26 

5.95 

0.78 

47.4 (116) 

III 

MS 

100 

25 

0.0160 

15.0 

0.1 



0.7 



1.26 

5.39 

1.06 

91.9 (261) 

IV 

KR 

100 

25 

0.0160 

5.5 / 1.8 

0.2 

/ 

0.1 

1.2 

/ 

0.4 

1.26 

5.16 

1.05 

111.7 (335) 

V 

SC 

100 

25 

0.0160 

5.8 / 2.0 

CO 

o 

/ 

0.1 

1.3 

/ 

0.4 

1.26 

5.37 

1.10 

89.9 (256) 

VI 

SA 

250 

62 

0.0116 

15.0 

0.1 



1.0 



1.71 

8.11 

1.40 

118.1 (288) 

VII 

TA 

250 

62 

0.0116 

15.0 

0.1 



1.0 



1.71 

8.11 

1.32 

184.6 (450) 

VIII 

MS 

250 

62 

0.0116 

15.0 

0.1 



0.6 



1.71 

7.03 

1.37 

171.8 (517) 

IX 

KR 

250 

62 

0.0116 

5.5 / 4.8 

0.2 

/ 

0.1 

0.9 

/ 

0.4 

1.71 

6.75 

1.66 

177.5 (568) 

X 

SC 

250 

62 

0.0116 

5.8 / 5.1 

0.2 

/ 

0.1 

1.0 

/ 

0.5 

1.71 

7.20 

1.62 

217.8 (634) 

XI 

SA 

500 

125 

0.0123 

15.0 

0.1 



1.0 



2.15 

10.19 

2.00 

205.1 (498) 

XII 

TA 

500 

125 

0.0123 

15.0 

0.1 



1.0 



2.15 

10.19 

1.91 

280.0 (680) 

XIII 

MS 

500 

125 

0.0123 

15.0 

0.1 



0.6 



2.15 

8.75 

1.98 

204.4 (625) 

XIV 

KR 

500 

125 

0.0123 

5.5 / 5.1 

0.2 

/ 

0.1 

0.9 

/ 

0.4 

2.15 

8.51 

2.02 

270.1 (859) 

XV 

SC 

500 

125 

0.0123 

5.8 / 5.4 

0.2 

/ 

0.1 

1.1 

/ 

0.5 

2.15 

9.05 

1.98 

262.4 (763) 

XVI 

SA 

750 

187 

0.0081 

15.0 

0.1 



1.0 



2.47 

11.68 

2.30 

329.5 (804) 

XVII 

TA 

750 

187 

0.0081 

15.0 

0.1 



1.0 



2.47 

11.68 

2.30 

387.5 (945) 

XVIII 

MS 

750 

187 

0.0081 

15.0 

0.1 



0.6 



2.47 

9.98 

2.32 

162.6 (503) 

XIX 

KR 

750 

187 

0.0081 

5.5 / 5.2 

0.2 

/ 

0.1 

0.9 

/ 

0.4 

2.47 

9.73 

2.27 

227.0 (729) 

XX 

SC 

750 

187 

0.0081 

5.8 / 12.9 

0.2 

/ 

0.1 

1.0 

/ 

0.6 

2.47 

10.30 

2.27 

369.8 (1061) 


* SA Salpeter IMF, TA Taff IMF, MS Miller & Scalo IMF, KR Kroupa IMF, SC Scab IMF. 

* Total number of stars (N s + 2 N b ). 
x Semi-major axis for PBs in pc. 

t In AIq. For KR and SC models (Binary / Single). 

* Initial virial radius in pc. 

* Initial tidal radius in pc. 

° Initial half-mass radius in pc. 

0 Disruption time in scaled units (in Mys). 


complex interplay between the IMF, mass loss from stel¬ 
lar evolution and primordial binaries. For poor clusters the 
massive stars in PBs dominate the cluster evolution. For 
rich clusters, the interplay between power-law IMFs, PBs 
and stellar evolution seems to increase the cluster life-time 
in relation to models with PBs but no stellar evolution. 
For intermediate population (N = 250, 500), power-law 
models with all realistic features disrupt earlier than their 
respective models without mass loss. However, there is no 
clear interpretation for this behaviour. The increased life¬ 
time in some models for N = 750 can be explained by the 
formation of temporarily stable multiple systems (in some 
cases, hierarchical triple systems) in the cluster remnant; 
for N = 500, a certain number of these systems are also 
formed in some models. These systems delay the disrup¬ 
tion of the evolved cluster because some of them are long- 
lived. Formation of such systems depends strongly on the 
fraction of PBs, the cluster membership and the param¬ 
eter 1Z. Larger N and PB fraction promote an increased 
probability of formation of such multiple systems. 


3. Characteristic quantities 

In this section we compare all of the different runs, con¬ 
centrating, as in previous papers, on two representative 
diagnostic parameters. The first one measures the degree 
of dynamical evolution of the entire cluster. This quan¬ 
tity is called the evolution modulus and is defined by (von 
Hoerner 1976) 

W = log(R h /R c ), (5) 

where Rh is the half-mass radius and R c is the core radius. 
The core radius defined by Casertano and Hut (1985) is 

4orb' = ^ Pi Ri /Si Pi . (6) 

Here, Ri is the distance from star i to the density center, 
defined as the density-weighted centroid of the system, 
the density pi is determined by the distance Reg to star 
V s sixth nearest neighbour: 

Pi = Meg/Reg , (7) 

where Meg is the total mass lying within distance Reg of 
star i (excluding the mass of star i itself), and the sum is 
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taken over all stars in the system. The second parameter 
is the escape rate dN/dt which describes the disruption 
rate of the system. 

3.1. Evolution Modulus 

This quantity has the advantage that it can be obtained 
directly from observable properties. Figs. show the 


Fig. 2. Evolution modulus for Salpeter’s models with N = 
100, 250, 500, 750. In all figures, time is given in units of the 
half-mass initial crossing time in the current system. 


Fig. 3. Evolution modulus for Taff’s models with N = 100, 
250, 500, 750. 

behaviour of W as a function of time for all the models. 
If we compare the present figures with those from Paper 
II we observe significant differences. 

For Salpeter IMF models, the behaviour of W is very 
similar for all N in Paper II, however the inclusion of PBs 
shows great differences, depending on N, for our present 
models. First, the evolution time-scales (in scaled units) 
are very different depending on N, accelerating the dy¬ 
namical evolution in poor clusters. For N = 100, the clus¬ 
ter is disrupted even having W > 0. For N = 250 the 


Fig. 4. Evolution modulus for Miller & Scalo’s models with N 
= 100, 250, 500,750. 


Fig. 5. Evolution modulus for Kroupa’s models with N = 100, 
250, 500, 750. 

evolution is also accelerated but for N = 500 the evolu¬ 
tive time-scales are almost the same as without PBs. The 
evolution is even slower for N = 750. As before, this trend 
can be explained by the formation of temporarily bound 
multiple systems (triple and quadruple). The comparison 
between our present results for W and those from Paper 
III also presents several differences. The evolution of W is 
slower for greater N but for smaller N it is nearly similar. 
The mean value of W is now significantly smaller because 
the halo is less extended than in the case of models with 
PBs but no stellar evolution. 

For Taff IMF models, the behaviour of W is roughly 
similar to that presented in Paper II except in the case of 
N = 100. In this case, the cluster evolution is accelerated 
in a very significant way, but the cluster is completely dis¬ 
rupted after reaching W = 0. The behaviour is different if 
comparing with Paper III because the core-halo structure 
seems to be stabilized by stellar mass loss as compared 
with models with PBs but no mass loss. 

Models with Miller & Scalo IMF show practically sim¬ 
ilar behaviour to those presented in Paper II. The only 
main differences appear for early stages of the cluster evo¬ 
lution where the mean value of W is greater for models 
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Fig. 6. Evolution modulus for Scalo’s models with N = 100, 
250, 500, 750. 

without PBs. This suggests that the haloes formed in mod¬ 
els with PBs and stellar evolution are less extended. As 
regards comparison with models without stellar evolution 
we also observe features suggesting that the interplay be¬ 
tween mass loss and PBs produces less energetic haloes 
than in the case of models with PBs but no stellar evolu¬ 
tion. 

For models with Kroupa IMF, we observe that the evo¬ 
lution is very different from that presented in Paper II. For 
large N the evolution is slowed down but for N = 100 the 
evolution is accelerated although less than for power-law 
models. The evolution of W is affected even at early stages 
of the cluster evolution, when the stellar mass loss is not 
yet dominant. If comparing with Paper III, the evolution 
of the cluster seems to be slightly accelerated. 

Finally, for Scalo IMF we observe that with regard to 
Paper II the evolution for models with N > 250 is signifi¬ 
cantly slowed down, but speeding up a bit for smaller N. 
As regards Paper III, the ratio core-halo (in size) seems to 
be more stable for models with mass loss and the overall 
evolution is slowed down. 

We observe that the inclusion of stellar evolution in 
models with PBs affects their evolution considerably but 
in a very uncertain way. However, it is clear that the 
changes are very IMF dependent. The interpretation of 
the results for our present models in relation to models 
with PBs but no stellar evolution is very unclear in com¬ 
parison with the conclusions we obtained with regard to 
models without PBs but with or without mass loss due to 
stellar evolution. 

3.2. The escape rate 

This is another important quantity for the cluster evo¬ 
lution. Figures 0-0 show the number of cluster stars as 
a function of time. For N = 100 the inclusion of stel¬ 
lar evolution in models with PBs distinguishes clearly be¬ 
tween models with power-law IMFs and realistic ones. In 
comparison with figures for the escape rate without stellar 


Fig. 7. Evolution of the cluster population with time for mod¬ 
els with N = 100. 


Fig. 8. Evolution of the cluster population with time for mod¬ 
els with N = 250. 


evolution, the disruption of power-law models is acceler¬ 
ated but in the case of realistic IMFs it is slowed down. 
Even for Scalo IMF we observe an increase in the dura¬ 
tion of the stage in which close encounters are dominant. 
The almost exponential behaviour that we observe in some 
models in Paper III disappears because most of the esca¬ 
pers are not due to distant encounters. As regards com¬ 
parison with Paper II the evolution of the escape rate is 
very different so we suggest that the escape mechanism 
induced by binaries with stellar evolution is different from 
the dominant one in clusters with mass loss due to stellar 
evolution but no PBs. 

The same trend is also observed for N = 250 but it is 
not so clear as in the case of IV = 100. Models with Taff 
and Miller & Scalo IMFs show a nearly similar behaviour 
as in Paper III. 

For N = 500 the behaviour of the escape rate is 
very different from that presented in Paper III. Models 
with massive stars (power-law IMFs) lose their distinctive 
features and accelerate their disruption rates (in scaled 
units). For the other IMFs the figures are nearly simi¬ 
lar to their respective ones without PBs. With regard to 
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Fig. 9. Evolution of the cluster population with time for mod¬ 
els with N = 500. 


Fig. 10. Evolution of the cluster population with time for 
models with N = 750. 

comparison with Paper III, the escape rate is reduced by 
including stellar evolution. 

For IV = 750 the inclusion of stellar evolution in models 
with PBs produces a significant change in the escape rate 
for power-law IMFs. Their escape rate show a slow-down 
with respect to those from Paper III. However, models 
with realistic IMFs, except for Scalo IMF, look very sim¬ 
ilar to those from Paper III. With regard to comparison 
with Paper II results, the escape rate is slowed down for 
models with Miller & Scalo, Kroupa and power-law IMFs 
but is accelerated for Scalo IMF. 

As we can see, the interpretation of the results is diffi¬ 
cult in most of the cases. Only for N = 100 can we do this 
in an easy way. Clearly the differences induced in these 
models are affected by the massive stars. The stellar evo¬ 
lution in models with small N and PBs is the dominant 
mechanism for the dynamical evolution of open clusters. 
For most of the models with small N we found in Pa¬ 
per III almost exponential decay in the escape rate which 
suggests an evaporative (by distant encounters) dominant 
escape mechanism, however now we find a change in the 
shape of the curve. The reason seems to be due to the halo 
extension. Models including stellar evolution seem to show 


less extended haloes than models without, so the number 
of stars available for leaving the system in a smooth way 
is smaller. This trend almost dissapears for models with 
N = 250 but it can still be observed preferentially at early 
stages of the cluster evolution. For larger N the explana¬ 
tion is not very clear, also we note less extended haloes but 
the results depend strongly on the IMF. In any case the 
dominant mechanism for the escape rate at early stages of 
the star cluster evolution seems to be by close encounters. 

4. Evolution of PB population 

This section is devoted to study the evolution of the PBs 
in our present models. We are mainly interested in com¬ 
paring with results from Paper III. In Figs. [i~T|4l~~i| we see 
the percentage of primordial binaries in the cluster as a 
function of time. 


Fig. 11. Evolution of the PB population as a function of time 
for models with N = 100. 


Fig. 12. Evolution of the PB population as a function of time 
for models with N = 250. 

From the figures, the first thing to note is the life time 
of the primordial binary population. For poorly populated 
models (N = 100), we obtain that the life time of the PBs 
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Fig. 13. Evolution of the PB population as a function of time 
for models with N = 500. 


Fig. 14. Evolution of the PB population as a function of time 
for models with N = 750. 

are significantly shorter than for Paper III models. More¬ 
over, the evolution of the surviving fraction of the PB 
population for models with power-law IMFs seems to be 
clearly different from that observed for models with real¬ 
istic IMFs. The diminution in the number of PBs is more 
rapid because the evolution of the cluster itself is quicker. 
The reason is due to the supernova events. In Salpeter 
and Taff models two supernova events have occurred; as a 
result the cluster disintegration and the escape or disrup¬ 
tion of PBs are accelerated. For N = 250, some differences 
are observed as regard models from Paper III. Salpeter 
models show a significant acceleration in the diminution 
rate of the PB fraction caused, as before, by the super¬ 
nova explosions (2 for these models). However, the Taff 
model shows a nearly similar behaviour as in Paper III; 
the same can be observed for realistic IMFs except for the 
Scalo model. For N = 500 the almost exponential decrease 
in the percentage of surviving binaries which appears in 
models from Paper III disappears for the majority of the 
present models. Initially the same linear diminution is ob¬ 
served, but later the diminution rate slows down by a very 
significant value. For some models in Paper III the per¬ 
centage of surviving PBs at the cluster mean life time is 


smaller than 30 % but now this value reaches almost 50 
% for power-law IMF models. Models with N = 750 also 
show differences. In Paper III, all the models (excluding 
the monocomponent one) have a nearly similar behaviour 
but now models with Salpeter, Taff and Scalo IMFs show 
almost the same behaviour but the others do not. For the 
former models the diminution of PBs is almost linear and 
smaller than for Miller & Scalo and Kroupa IMFs. From 
Table [| we see that models for Salpeter, Taff and Scalo 
IMFs have more massive stars so this may be the reason 
for this behaviour. In average, our models show that PBs 
are retained preferentially; this result has also been found 
by Aarseth (1996d). 

Our present considerations suggest that the member¬ 
ship is a main factor for retaining PBs; models with in¬ 
creasing N show smaller diminution rates for the percent¬ 
age of surviving PBs. Also, the IMF has a main role in 
the evolution of PBs; for poorly populated clusters mas¬ 
sive stars in binaries control the evolution of the entire 
cluster. For small open clusters the presence of massive 
stars in binaries can accelerate their disruption but for 
more populated star clusters the effect is to the contrary. 


Table 3. Binary fraction in the core and in the whole cluster 


MODEL f° 

ft 

fc 

ft 

fc 

ff 

I 

0.31 

0.33 

0.00 

0.35 

0.00 

0.13 

II 

0.33 

0.33 

0.00 

0.28 

0.00 

0.57 

III 

0.63 

0.33 

0.80 

0.28 

0.36 

0.57 

IV 

0.45 

0.33 

0.33 

0.27 

0.18 

0.57 

V 

0.47 

0.33 

0.27 

0.23 

0.09 

0.22 

VI 

0.33 

0.33 

0.32 

0.31 

0.25 

0.50 

VII 

0.35 

0.33 

0.00 

0.31 

0.00 

0.83 

VIII 

0.30 

0.33 

0.33 

0.31 

0.27 

0.57 

IX 

0.33 

0.33 

0.17 

0.29 

0.09 

0.57 

X 

0.36 

0.33 

0.36 

0.39 

0.33 

0.71 

XI 

0.36 

0.33 

0.75 

0.34 

0.19 

0.57 

XII 

0.37 

0.33 

0.42 

0.31 

0.09 

0.38 

XIII 

0.34 

0.33 

0.50 

0.32 

0.09 

0.83 

XIV 

0.36 

0.33 

0.33 

0.28 

0.08 

0.63 

XV 

0.27 

0.33 

0.33 

0.28 

0.25 

0.33 

XVI 

0.34 

0.33 

0.00 

0.35 

0.00 

0.63 

XVII 

0.33 

0.33 

0.00 

0.28 

0.00 

0.80 

XVIII 

0.32 

0.33 

0.47 

0.32 

0.33 

0.50 

XIX 

0.37 

0.33 

0.55 

0.25 

0.13 

0.22 

XX 

0.36 

0.33 

0.55 

0.30 

0.29 

0.56 


An interesting parameter to study is the binary frac¬ 
tion, /, both in the core and the whole cluster. Table || 
gives the core and total binary fractions for all the models 
at three selected times in the cluster life: at time equal to 
zero (/,?, / t °), at the time in which the stellar population 
is a half of the initial one (fjf, ft ) and at the end of the 
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simulation (/®, / t e ). In this Table we only consider the pri¬ 
mordial binaries because they are the only hard ones; in 
our present models binaries which are formed dynamically 
are all soft and short lived. Some quick conclusions arise 
from the present values of the binary fraction. First, for 
power-law models there is a strong tendency for binaries 
to underpopulate the core on a short time scale. Except for 
N = 500 all the power-law models show a smaller binary 
fraction in the core at the half cluster life. The explana¬ 
tion is due to the supernova events located in the core. 
They dominate the evolution, driving the binaries outside 
the core (values for the total binary fraction are nearly 
similar for all the models at the cluster half life). A typi¬ 
cal behaviour of the total binary fraction can be observed 
in Fig. [l^. In all the models an initial decay is observed 
but the total binary fraction increases in most of the cases 
when the cluster population has decreased significantly. It 
suggests a preferential escape of the single stars in clus¬ 
ters with a fraction of PBs. However, the evolution of the 
core binary fraction with time is extremely irregular. Its 
behaviour is extremely chaotic although almost in all the 
models an initial increase is observed as in Aarseth (1996b, 
c). This effect is greater for rich models and realistic IMFs. 
Second, the value of the total binary fraction towards the 
end of our simulations for N < 250 is 0.57 in many of 
the models; it is the value currently accepted for the bi¬ 
nary fraction in the solar neighbourhood (Duquennoy & 
Mayor 1991). This value appears preferentially for models 
with Miller & Scalo and Kroupa IMFs so it is possible to 
suggest that there probably exists a link between the size 
of the cluster in which stars are born, their IMF and the 
binary fraction in a certain region of space. 


Fig. 15. Evolution of the total binary fraction as a function 
of time in years for models with N = 750. In all the models an 
initial decay is observed but it increases in most of the models 
when the cluster population has reached values smaller than 
one third of the initial. 


5. Velocity distribution of the escaping stars 

In this section we analyze the velocity distribution of the 
escaping stars, both singles and binaries. The main em¬ 
phasis is on comparing present results with previous ones 
for models without PBs and mass loss, with mass loss but 
no PBs and with PBs but no mass loss. Most of the results 
that appear in this section were not discussed previously 
in the papers of this series. 

Stars are considered as escapers and they are removed 
from the calculations when they reach a distance from 
the cluster centre greater than twice the tidal radius. At 
that moment, we note the velocity of the star (if single) 
or of the centre of mass (if binary). The escape bound¬ 
ary depends on the Galactic parameters so the obtained 
velocities are higher that those obtained from the theo¬ 
retical expression (2 < v >); many of our velocities have 
excess escape energy. The values of the mean escape ve¬ 
locity, the standard deviation in the mean, the minimum 
escape velocity and the maximum escape velocity can be 
seen in Table [J. We observe two clear tendencies: first, 


Table 4. Escape velocity for the models* 


MODEL 

< C Vesc ^ 

0~m 

min 

v esc 

max 
”esc 

I 

0.935 

0.070 

0.261 

4.450 

II 

0.866 

0.062 

0.413 

4.426 

III 

0.822 

0.056 

0.309 

3.511 

IV 

0.699 

0.044 

0.000 

2.413 

V 

0.710 

0.052 

0.223 

3.137 

VI 

1.190 

0.065 

0.251 

7.551 

VII 

1.270 

0.051 

0.434 

5.992 

VIII 

0.996 

0.044 

0.000 

5.768 

IX 

0.858 

0.030 

0.270 

3.994 

X 

1.049 

0.049 

0.210 

6.915 

XII 

1.419 

0.050 

0.000 

11.261 

XIII 

1.456 

0.057 

0.000 

18.427 

XIV 

1.140 

0.031 

0.245 

5.050 

XV 

1.051 

0.023 

0.241 

4.726 

XVI 

1.211 

0.036 

0.264 

7.819 

XVII 

1.558 

0.030 

0.345 

6.960 

XVIII 

1.544 

0.036 

0.345 

8.687 

XIX 

1.141 

0.024 

0.267 

6.386 

XX 

1.158 

0.026 

0.000 

7.988 

XXI 

1.246 

0.030 

0.000 

11.080 


* All the data in km/s. 

the mean escape velocity increases with N and second, 
the dispersion of velocities decreases when N increases. 
This is to some extent connected with the physical scaling 
in the models. Moreover, the greater mean escape veloc¬ 
ities are always for models with Salpeter and Taff IMFs, 
so it depends on the IMF. This is to be expected because 
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the masses of these models are greater than for the other 
IMFs. As we can see, when N increases the probability of 
high velocity escapers grows; we observe that in all cases 
the maximum escape velocity is considerably higher than 
the mean escape velocity. The origin of these high veloc¬ 
ity escapers is connected with binary-binary interactions 
so the high escape velocities must also be connected with 
binaries of a certain size. 

As regards comparison with results from previous 
models, we observe significant differences. Considering re¬ 
sults from models of Paper III, the mean velocity of the 
escapers is reduced in most of the present models. This 
tendency is clear for clusters with simple power-law IMFs. 
However, realistic IMFs show a more complex behaviour. 
For N = 100, the same tendency as for old IMFs is ob¬ 
served (except for Kroupa IMF). For N = 200, the mean 
velocity is larger for the current models. For N = 500, this 
is only observed for Miller & Scalo IMF. For N = 750, the 
same tendency is observed for all the models. Concerning 
the theoretical implications of these experimental results, 
it seems that stellar evolution of massive stars is the main 
process which explains the observed differences. Mass loss 
in binaries promotes the disruption of a certain number 
of such systems for models with simple power-law IMFs 
so the global results suggest the number of initial binaries 
would have been smaller. On the other hand, present val¬ 
ues for the mean velocity of the escapers are larger than 
those of Paper II, so the presence of PBs increase the ve¬ 
locity of the escapers (due to close encounters between 
binaries and singles) for models with mass loss. 

As regards the relation of velocities of escaping stars 
with their masses, we observe a bimodal distribution. Most 
of the escaping particles have velocities not much larger 
than the mean but a few per cent has velocities greater 
than three times the mean. This suggests two processes for 
the escape of cluster stars. Low velocity escapers are gen¬ 
erated by evaporation, i.e. gradual increase of the velocity 
because of distant encounters and high velocity escapers 
are produced by close encounters between singles or, more 
frequently in models with PBs, between singles and bina¬ 
ries or binary-binary encounters. 


6. Mass loss in binaries 

In order to obtain the change of mass as a function of time, 
NBODY5 includes the fast fitting functions by Eggleton 
et al. (1989) for Population I stars. For convenience, mass 
loss is implemented at discrete intervals when the accu¬ 
mulated contribution reaches 1 percent. The actual mass 
loss is assumed to be instantaneous, with no further effect 
on the cluster members or binary companions. Approxi¬ 
mate energy conservation can be achieved by performing 
appropriate corrections to the total potential energy. In 
the version of NBODY5 used, the binding energy is cor¬ 
rected by the term A M/r where A Ad is the mass loss 
and r is the binary separation. For close binaries the mass 


loss is always at apocentre, hence the effect will be small¬ 
est. When the mass loss is large, the binary is usually 
disrupted. 

The loss of mass from a binary system is a very 
complicated dynamical problem. Stellar mass loss causes 
changes in the orbital elements of the binary. The problem 
has been mainly studied by Hadjidemetriou (1963, 1966, 
1968). 

The orbital period can increase or decrease secularly 
depending upon the mass-flow conditions. The simplest 
case is where mass is lost isotropically from the system. 
This situation is assumed in our stellar mass loss events. 
In our models, the vast majority of massive binaries are 
disrupted during mass loss events, so the preferential bi¬ 
nary survivors are low mass binaries in which the mass loss 
rate is sufficiently slow to permit a quasi-stable evolution. 

By Kepler’s third law, 


4 7r 2 a 3 
P 2 


G M 


( 8 ) 


where P is the binary period. It gives, for a constant semi¬ 
major axis a, the following relation between the change in 
period A P for a mass loss A M: 


AP _ AM 
~P~ ~ ~ 2 M ' 


(9) 


From Eq. (|tj) an abrupt change in the period can be 
achieved by one component losing mass in an eruptive out¬ 
burst, with the lost material ejected at a high speed. From 
our simulations, it is found that the evolution of surviving 
low mass PBs (not exchanged) with the two components 
of the same mass can be described almost exactly by Eq. 
(J)j) with small deviations due to the perturbers; i.e. the 
semi-major axis remains almost constant during most of 
the time when the mass loss rate is not too high. The 
evolution of massive binaries is very complicated because 
of mass loss. The majority of these binaries are destroyed 
during the first few crossing times. Their products, mainly 
white dwarfs, appear sometimes as members of another bi¬ 
nary. These binaries containing collapsed objects have a 
mass ratio of nearly one. 

The above formulae can be used for estimating the 
changes in the period for an unperturbed binary. The phe¬ 
nomena is more complex if we consider that the binary has 
several (in some case many) perturbers. In this case mass 
loss events near apocentre can easily disrupt the binary 
because the perturbing force is of the same order as the 
force between binary components. This depends mainly 
on the binary semi-major axis; binaries with small semi¬ 
major axis are like single stars from a dynamical point of 
view, so the mass loss events are less destructive for them. 
On the other hand stellar mass loss in clusters with PBs 
can favour the survival of certain kinds of binaries; bina¬ 
ries with both low-mass primary and secondary do not 
suffer a significant mass loss before their possible escape 
from the cluster. There is observational evidence (Verbunt 
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et al. 1994) that favour the above hypothesis; most X-ray 
sources in the old open cluster M67 are RS CVn binaries. 
These are active binary systems in which the primary is 
a star with a spectral type F or G in the main sequence 
or subgiant (luminosity class V or IV) and the secondary 
is a bit cooler and usually more massive and evolved with 
spectral type K and class IV. Their orbital periods are 
usually in the range 1 14 days. In models for N = 750 
we have found several binaries of this spectral type in the 
remnant, typically with a period of a few hundred years 
although in other models (not detailed here) with N = 
1500 and 500 primordial binaries the typical value for the 
period is about 8 days. In any case, the mean value for 
the periods of the surviving binaries depends on the ini¬ 
tial period distribution. 

7. Evolution of the stellar content 

As pointed out above all our models start with the stars 
on the ZAMS. This assumption implies some questionable 
hypotheses. First, it assumes that the stellar formation 
in the cluster took place by a single burst. On the other 
hand, it means there is no preferential stellar mass for 
the beginning of the star formation. The second of these 
premises is arguable in the light of some observational 
discoveries (Herbig 1962; Iben & Talbot 1966; Cohen & 
Kuhi 1979; Adams et al. 1983; Strom 1985). They suggest 
that low-mass stars form first and over a long time scale. 
However, Strom (1985) cautioned that these conclusions 
rested heavily on the theoretical pre-main-sequence tracks. 
As pointed out in the excellent review by Zinnecker et 
al. (1993), low-mass stars acquire their final mass first, 
but will then take a long time to reach the ZAMS; high- 
mass stars take longer to accumulate their mass but they 
evolve rapidly onto the ZAMS. Indeed, the data available 
are consistent with simultaneous formation of stars of all 
masses (Stahler 1985; Schroeder & Cornins 1988). In spite 
of these questionable initial hypotheses, the evolution of 
the stellar content of our models is analyzed in the present 
section. 

We want to study specially the destiny of the collapsed 
objects in our models. As we see from Table ||. some of our 
models (preferentially with power-law IMFs) have stars 
with enough mass to finish their nuclear life in a super¬ 
nova event. From our current knowledge of stellar evo¬ 
lution, after such an episode the final product may be a 
neutron star (called pulsar if its radio emission can be de¬ 
tected from the Earth) or a black hole. The last kind of 
collapsed object has not been introduced in the stellar evo¬ 
lution routines used in our present models. As regards the 
formation of neutron stars, it is considered in our models. 
There are a significant number of neutron stars (pulsars) 
detected in globular clusters, but at present there is no 
evidence for these collapsed objects in open clusters. In 
our models there are no surviving neutron stars. They es¬ 
cape from the cluster shortly after their formation. It is 


because after the formation of a pulsar, it suffers a kick ve¬ 
locity due to a strong non-symmetric mass loss. Although 
models with greater N can retain a population of neutron 
stars, our present models are not able to keep the neutron 
stars. The time-scale for leaving the cluster for a newly- 
born neutron star in our models is a few kyr, so it can 
be considered as an explanation for the null population 
of neutron stars observed in open clusters. However, it is 
to be expected from models with greater N (N = 1500 in 
models not presented here) that there is a higher (but low) 
chance of detecting a neutron star in a highly populated 
open cluster (maybe M67 is a good candidate to achieve 
this). 

On the contrary, white dwarfs remain in the cluster 
for many crossing times in most of the models. We have 
even found a few models in which there are one or two 
white dwarfs among the members of the cluster remnant. 
However, there are some differences between the life times 
for white dwarfs in our models depending both on the IMF 
and cluster membership. 

8. Open Cluster Remnants 

In this section we consider the composition of the final 
cluster remnant in our models. In Paper III, we noted 
that the final remnants had a common feature, a high bi¬ 
nary richness. This is also observed in the present models 
but the fraction of purely primordial binaries remaining in 
the final object resulting from the cluster disruption has 
decreased. Also we observe that the binaries with two com¬ 
ponents of the same mass are the preferential survivors. 
In the majority of the models the number of remaining 
binaries is three or four (we stop the simulations when 
the cluster population is N = 10). From the results of 
models for Papers I and II, we also find a certain number 
of binaries in the cluster remnant (usually 1 or 2). These 
binaries are formed dynamically, not primordials so one 
of them is hard and the others (if any) are soft. In our 
present models all of the surviving binaries are hard. 

The final binary population shows distinctive char¬ 
acteristics depending on the cluster richness. For mod¬ 
els with N < 250 the surviving binaries do not have a 
preferential ratio between the masses of the binary com¬ 
ponents. We find binaries with both components being 
low-mass stars (A4 < O.8M 0 ) or both massive or one of 
them massive and the other not. However, for rich models 
the surviving binaries are always with almost equal-mass 
components. In a few cases one component is a giant and 
the other is a low-mass star with a spectral type typically 
in the range M0-M5. When the two components are on 
the ZAMS both the primary and secondary have spectral 
type in the range K5-M5. In general, the stellar content 
of the open cluster remnants is mostly ZAMS stars and 
in a few models there appears a white dwarf or a red gi¬ 
ant, especially for rich clusters. In most of the models the 
binary fraction in the remnant is significantly larger than 
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the initial one; Aarseth (1996d) has also found this effect. 
In a model with 10,500 stars and 500 binaries, the remnant 
contains 10 binaries and 39 stars. 

Although binaries in cluster remnants are primordial, 

i.e. at least one of their components was in a PB at the be¬ 
ginning of the simulation, not all the surviving binaries are 
purely primordial. For rich clusters we find a certain frac¬ 
tion of exchanged binaries (1 or 2 per model). For small 
clusters it is rare to find an exchanged binary because the 
cluster disrupts faster so the probability of an exchange 
during the cluster life is reduced. 

The explanation of the above results rests on the time- 
scales for stellar evolution. For clusters with small N, the 
cluster disrupts before massive stars have started to leave 
the ZAMS, so we can find binaries with massive primaries 
in the cluster remnant. However for rich clusters, the dis¬ 
ruption time-scale is greater than the characteristic time- 
scale for significant mass loss from massive stars, so bina¬ 
ries with massive primaries could be disrupted or ejected 
during mass loss events. The observational test of our re¬ 
sults is very difficult in the case of rich clusters remnants. 
The detection of a small population (a few tens) of faint 
stars most of them less luminous than our Sun is a big 
challenge for the biggest optical telescopes. Moreover, the 
detection of the binaries in these rich cluster remnants is 
even more difficult. However, there is some light from this 
dark picture because stars with such spectral types are 
expected to be strong radio and/or X-rays sources which 
can be detected with synthesis techniques or by instru¬ 
ments on satellites. The search of open cluster remnants 
has been considered by Loden (1977, 1979, 1980), who 
analyzed very loose or star-poor clusterings in a large sur¬ 
vey of the Southern Milky Way. Thousands of these ob¬ 
jects were found and classified into 4 sets. One of them is 
formed by extremely small and star-poor clusters which he 
suggested as possibly cluster remnants. The frequency of 
these objects in his sample is about 20 % but he considers 
that it can be significantly greater. 

9. Conclusions 

The main conclusions from this work can be summarized 
as follows: 

1. - The inclusion of stellar evolution in cluster mod¬ 
els with a fraction of PBs affects the overall dynamical 
evolution of the cluster in a uncertain way which depends 
strongly on the cluster richness and the IMF. 

2. - The stellar evolution in clusters with small popula¬ 
tion and power-law IMF accelerates their disruption in a 
very significant way. 

3. - The escape velocity increases with the cluster rich¬ 
ness but the dispersion decreases when N increases. Some¬ 
times, a star can leave the cluster with high velocity. 

4. - The final cluster remnant is very rich in binaries, 
frequently purely primordial but its composition depends 
strongly on the initial cluster population because of the 


interplay between the time-scales for cluster disruption 
and stellar mass loss. Binaries in remnants of poor clusters 
do not have any special feature in their components but 
binaries in rich cluster remnants have usually almost the 
same mass for their components and have late spectral 
types. Collapsed objects are almost always absent from 
open cluster remnants. 
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